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Abstract. In this paper we propose an improved method 
for calculating Henon's stability parameter, which is based 
on the differential of the Poincare map using the first vari- 
ational equation. We show that this method is very accu- 
rate and give some examples where it gives correct results, 
while the previous method could not cope. 
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1. Introduction 

Over 30 years ago, in a now seminal paper, Henon 
(1965) introduced the stability parameter a, which distin- 
guished whether a given periodic orbit is stable or not. 
This distinction is essential since the properties of the 
two types of orbits differ considerably. Indeed stable peri- 
odic orbits trap regular orbits around them, while unstable 
ones trigger chaos. The use of Henon's stability parame- 
ter is, however, not limited to that. By allowing us to find 
precisely the value of the energy for which a given family 
changes from stable to unstable, or vice-versa, it allows us 
to find the point from which new families bifurcate, since 
for a periodic orbit a transition from stability to instabil- 
ity results in a bifurcation of a new stable family, while 
a transition from instability to stability introduces a new 
unstable family. 

To calculate stability, Henon (1965) approached the 
differential of the Poincare map by finite differences. This 



technique has been so far widely followed in galactic dy- 
namics (e.g. Contopoulos & Gr6sb0l 1989 and references 
therein). Nevertheless it suffers from a number of disad- 
vantages. We have found it quite adequate in regular re- 
gions, but found it could not cope with difficult chaotic 
regions. For this reason we present here an alternative ap- 
proach, based on the differential of the Poincare map using 
the first variational equation. We will hereafter refer to it 
as variational equation method. The variational equation 
technique is used in other domains that need very accu- 
rate results, like the study of the three-body problem, the 
Stormer problem or a few others problems in celestial me- 
chanics (e.g. Deprit & Price 1965; Markellos 1974; Markel- 
los & Zagouras 1977), as well as to obtain the Lyapunov 
exponents (Benettin et al., 1976, 1980; Contopoulos et al., 
1978; Udry & Pfenniger 1988). In this paper the varia- 
tional equation will allow us to obtain the differential of 
the Poincare map exactly. In section ?? we give a very 
brief mathematical justification of this method, which is 
well known from other applications, and then we apply it 
to the particular case of a two-dimensional system which 
is stationary in a rotating frame of reference, a case of- 
ten encountered in galactic dynamics problems. Section 
?? gives an example demonstrating the accuracy of the 
method we propose and section ?? another example for 
which this new method gives accurate results, while finite 
differences can not cope. We conclude in section ??. 



2. Method 
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2.1. Notation 

Let us consider an autonomous dynamical system, ex- 
pressed by the ordinary differential equation (hereafter 
ODE). X = F{x) with x e ST and where x = ^. By 
a solution of this ODE we will mean a map (p : U C 
MxFT — > FT such that (i>{t, x) = F{(j){t, x)). 

Using the precepts of Poincare (1892) wc can reduce 
the study of a continuous time system (ODE) to the study 
of an associated discrete time system (map) called the 
Poincare map, P : V C T, — »• E; x i — > x = (j){T{x),x), 
where S is a hyper-surface perpendicular to the flow F, 
which we will hereafter refer to as section, V is an open 
set in E, and t{x) is the time necessary for the point x 
to return for the first time to the section with the same 
sense of transversal of E. With this technique the problem 
of calculating the stability of a periodic solution of a ODE 
is reduced to the problem of calculating the stability of a 
fixed point of a map. For this we must check if solutions 
starting close to a fixed point at a given time remain close 
to it for all later times (Lyapunov stability). We thus com- 
pute the Taylor series expansion and we study the linear 
term of such an expansion. 

Let xq be a fixed point of P and x = xo + Axq a point 
in its neighbourhood. The Taylor expansion of P{x) is 

P{x) =P{xo)+DP{xo)Axo + 0{A^) (1) 

If we denote P{x) = xq + Axi, we have the linear 
relation Axi = DP{xo)Axo. The eigenvalues of DP will 
determine the stability. In two dimensions, if P is an area 
preserving map, the periodic orbit (j){t, xq) is stable if |q;| < 
1, where a is the stability parameter introduced by Henon 
(1965) and defined as a = |(aii +022) where DP = (a^). 

2.2. Calculation of DP 

Henon (1965) approximated the elements of the Jacobi 
matrix DP{x) using finite differences, i.e. 

Z,PW.,^ ^fe + ^]-^fe' +0(A) (2) 

As will be shown below, the above approximation gives 
sufficient accuracy in regular regions, but not in regions 
dominated by chaotic dynamics. 
We can calculate DP exactly, as: 

DP{x) = ^{t{x), x)Dt{x) + D(I){t{x),x) (3) 



The matrix D(j) can be obtained as a solution of the first 
variational equation with initial condition D(p{0, x) = Id. 

j^D^it, x) = DF{<j>{t, x))D4>{t, x) (4) 

To compute Dt{x) we use the fact that we are on the 
hyper-surface. Defining the hyper-surface E = g(x) = 
and differentiating we obtain after some algebra: 

where xi = (f){T{x), x) and where the symbol ( , ) rep- 
resents the dot product. The denominator is different from 
zero since the hyper-surface E is perpendicular to the flow 
F. 

2.2.1. Particular case 

We will now apply the above general technique to a two 
dimensional system which is stationary in a rotating frame 
of reference. Let us consider an autonomous dynamical 
system, expressed by the following ODE 

X = + 2^6^ + ^Ix \ 
y = -^y- 2flbX + flly / 

where ^{x, y) is the potential, Ob is the pattern speed of 
the coordinate system in which the dynamical system is 
stationary, and and denote the partial derivatives 
of the potential with respect to x and y respectively. 

Writing the above second order ODE as a system of 
four first order ODE, and using the momenta X = x — Slb 
and F = y -|- Oft we get 

X = X + Q.by 1 
y^Y-ntx ( 

Y = -%- n^x J 

This can also be expressed in vectorial form i = ^(2:) 
with z = {x,X,y,Y) and since below we will need the 
individual components of F, we will denote them as Fj 
with j = 1 to 4. Hereafter, numbers as subindex indicate 
a component. 

Since on the Poincare section y = and since the en- 
ergy in a rotating frame, Ej{x,X,y,Y), is an integral of 
motion, we can express F as a function of x and X, i.e. 

Y = Y{x,X). This restricts the Poincare map to the two- 
dimensional space {x,X). since, however, the system of 
first order ODEs has four equations, the variational equa- 
tion and the method to calculate the differential of the 
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Poincare map hold in a four- dimensional space. We must 
therefore express the Poincare map in four dimensions in 
order to calculate the differential and finally project in 
the two-dimensional space {x, X) as shown schematically 
below 



ix,X,0,Y) 
ix,X) 



{x,X,0,Y) 
{x,X) 



> =>P = 7ro^oi 



(8) 



Where i includes {x,X) in using the section equa- 
tion and the integral of motion, ^ is the Poincare map 
in IR'^ and tt projects the two first coordinates in M^. We 
can therefore write the differential of the Poincare map as 
DP = Dtt o o Di. 

Let us now proceed as in the general case, taking the 
derivative of ^'(2;) 



D^{z) = ^{t{z), z)Dt{z) + D(P{t{z), z) 



(9) 



Here 4>{t{z),z) = F{(P{t{z), z)), D(I){t{z), z) is the solu- 
tion of the first variational equation with initial condi- 
tion D(p{0,z) = Id. To obtain Dt{z) we use the section 
equation y = 0, which in iR* is equivalent to ^'3(2) = 0. 
We differentiate it and we obtain that the third vector 
of the matrix equation (??) should be equal to zero, i.e. 
D'^3{z) = o, from which it follows that 

DcI>3{t{z),z) 



Dt{z) 



Fs{ct>ir{z),z)) 
Now we can express DP in matrix form 



(10) 



DP-- 



(11) 



ay*i \ / 1 \ 
1 \ 5^*2 5x*2 5y*2 1 
OlOojo 00 

\d^^4 dx^i dy^i J\Y^ Yx J 

Where Y is, as discussed above, considered as a function 

of X and X . We sec that in order to obtain DP we need to 
apply D"^ only to Di, that is, we need to apply D(j) only to 
Di. Then, instead of solving the variational equation with 
initial condition D(f){0,z) = Id, which gives us D(l){t,z), 
we solve it applied to the vectors in Di, which means solv- 
ing i} = DF{y) with initial condition w° = (1,0, 0,Kc)'^ 
and w° = (0, \,0,Yx)^ . We will now show how it can 
be computed in practice. 

We integrate simultaneously the orbit and the variational 
equation for the two vectors in Di. 



x = X + ClbV 
X = + QbY 
y = Y- ribx 
Y = -%- 



with (f){0, zq) 



f xo\ 




(12) 





1 






V3 




\V4J 


V 



fib 

-^xy fife 



1 



-Qb 
-^yx — ^6 —^yy 



with V = v° and with v = w'^ 

Finally with those vectors we write DP as 



\ 










V3 


J 





(13) 



DP = 



Fi -o\ 



-^2 _o -0 ^2 -0 



(14) 



Remember that vj indicates the first return to the section 
- in the same sense of traversal - after the point Vj . This 
method involves the integration of a system of 12 rather 
than 4 equations, because we don't move only the point, 
but also two vectors. This of course takes more compu- 
tational time, but it gives us an accurate value for DP, 
since the error in DP is of the same order as that of the 
section points. 

3. Application to an axisymmetric potential 

As a first example let us take the axisymmetric logarithmic 
potential 



(15) 



where Vq and Rc are constants taken, in our example, to 
be equal to 1. and 0.1 respectively (Binney & Tremaine 

1987). The orbits of the main, circular family present no 
difficulty, so both finite differences and the variational 
equation method give roughly the same results. We can, 
however, compare the accuracy of the two methods by cal- 
culating the determinant of DP, which we will hereafter 
refer to as D. Its elements for finite differences are given 
by eq. (??), while for the variational equation method by 
eq. (??). 

A perfectly accurate calculation would of course give 
D = 1, and for less accurate calculations the determinant 
will deviate more from unity than for more accurate ones. 
We calculated finite differences for Ax = 10^'', 10"'*, 10~^ 
and 10~^ and we found that 10~^ gives the most accurate 
results. They are compared with those of the variational 
equation method in Fig. ??, which shows for both the 
value of log|Z) — 1|. The difference in accuracy between 
the two methods is striking! The finite difference method 
gives an accuracy between 10~^ and 10~^, while the accu- 
racy of the variational equation method is bound only by 



4 



C. Barbera^'^ & E. Athanassoula^: On the calculation of the linear stability parameter of periodic orbits 



the accuracy of the orbit calculation. It is this limiting ac- 
curacy that results also in the "quantisation" of the result- 
ing values. Thus the variational equation method gives, in 
most cases, an accuracy of at least 10^° better than that 
obtained with finite differences. 

The improvement in accuracy depends on the case consid- 
ered and is more important in more chaotic regions. 



Using the variational equation method we were able 
to show that this family is stable for as long as we could 
follow it. On the other hand finite differences with Aa; = 
10"^ find that the orbits are unstable if Ej > -126283, 
while for Ax = lO"-*, Ax = 10'^ and Ax = 10~^ the 
values of the Jacobi constant for which the family changes 
from stable to unstable are -126150, -126000 and -126500 
respectively. Thus finite differences give results that are 
qualitatively different from those of the variational equa- 
tion method and that depend heavily on the adopted value 
of Aa;. With the help of Poincare sections we confirmed 
that indeed this family is stable and therefore that it is the 
variational equation method that gives the correct result. 
We next repeated the calculation again using finite differ- 
ences but this time keeping only the points for which the 
value of the determinant D is sufficiently close to unity and 
noted that in this way the erroneous values disappeared. 
Although we are thus able to remove erroneous values, we 
are not able to find the correct value for the stability pa- 
rameter, unless we use the variational equation method. 



Fig. 1. Value of the log \D — 1\ as a function of the Jacobi en- 
ergy Ej for the axisymmetric logarithmic potential. The values 
obtained with finite differences axe given by an open circle and 
those obtained with the variational equation method by a star. 



4. Potential of a barred galaxy 

In the previous section we discussed an example where 
the variational equation method gives quantitatively more 
accurate results. In other cases, however, the differences 
between the two methods can be even qualitative. 
For our second example we will use the model 1 of 
Athanassoula (1992, hereafter A92). For most principal 
families the results of the two methods are in rough agree- 
ment, although the accuracy of the variational equation 
method is always better by at least as much as what we 
saw in the previous example. However, for families whose 
characteristic curve approaches the curve of zero velocity 
asymptotically, the differences can be much more impor- 
tant and there can be disagreement even as to whether 
a given orbit is stable or unstable. As an example let us 
take a Lagrangian family, the second one from the right in 
the lower panel of Fig. 2 of A92. We calculated its stabil- 
ity using finite differences as well as using the variational 
equation method and compare the results in Fig. ??. 




-1.266x10^ -1.264x10^ -1.262x10^ -1.26x10^ -1.25BxlO^ 

Fig. 2. Stability parameter a as a function of the Jacobi energy 
Ej for a Lagrangian family calculated in the model 1 of A92. 
The results obtained with finite differences and a Aa; of 10~^, 
lO-"*, 10"^ and lO"*^ are shown respectively by open circles, 
filled squares, open diamonds and stars. Results obtained with 
the variational equation method arc shown with a solid line 
and filled triangles. In order to make the figure clearer wc have 
plotted symbols only at one out of every five calculated points. 
For finite differences we have only plotted points for which 
|D — 1| < .5. By applying a stricter criterion we can eliminate 
successively more points, but in this way we get no information 
on the stability of the higher energy orbits. 
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5. Conclusions 

In essence the difference between the two techniques is 
that finite differences involve the calculation of liniAK^o 
for every element in the Jacobian matrix of the Poincare 
map. Numerically it is not possible to calculate the 
limAs->o, so we have to calculate this expression for dif- 
ferent Ax in order to find the optimum one that gives 
the best approach to the limit. It is, however, well known 
that calculating the derivative from finite differences can 
be unstable, and that it may not be possible to obtain a 
good approximation of limAa-^o- Thus this method may 
need several trials (i.e. more computing time), an a poste- 
riori check of the results, and, in certain cases, it may not 
be possible to find an appropriate value of Ax. 

The new method here calculates the Jacobian matrix 
accurately and doesn't depend on any Ax. It may involve 
some extra computing time because it involves the inte- 
gration of more differential equations, but this is, more 
often than not, compensated by the fact that it needs 
only one trial. Its assets are that it involves no tunable 
parameter, needs no a posteriori checking of the results 
and its accuracy is only limited by the method used to 
find the periodic orbits. In regular regions the two meth- 
ods give similar results, although the variational equation 
method is always more accurate. In chaotic or "difficult" 
regions, however, e.g. regions involving small stable islands 
in a chaotic sea, using the variational equation method 
may prove essential for obtaining the correct value of the 
stability parameter. Furthermore, a study of such regions 
may be very time consuming and the fact that with the 
variational equation method one can find precisely the bi- 
furcation points of new families can be a big help. We 
thus recommend it for all calculations such as the ones 
presented here. 
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